function [G1,G2,result,result_LS_code_base,result_LS_code_rover]=calc_preprocessing(RNX_header_base, RNX_data_base, RNX_header_rover, RNX_data_rover, NAV_header, NAV_data, X0_base, epoch, options)
%% function [G1,G2,result,result_LS_code_base,result_LS_code_rover]=calc_preprocessing(RNX_header_base, RNX_data_base, RNX_header_rover, RNX_data_rover, NAV_header, NAV_data, X0_base, epoch, options)
%% Data organization and preparation before spp, DGNSS or phase processing
%%
%% Beilin Jacques - ENSG/DPTS - 2011-12-18
%% Clement Fontaine - 2013-12-17
%%
%% Input
%%	- RNX_header_base, RNX_data_base : obs rinex data generated by function load_rinex_o
%%       If spp : set RNX_header_base and RNX_data_base to an empty structure
%%	- RNX_header_rover, RNX_data_rover : obs rinex data generated by function load_rinex_o
%%  - NAV_header,NAV_data : nav data (rinex_n or sp3) generated by load_rinex_n (for broadcasted ephemeris) or load_sp3 (for precise orbits)
%% 		- broadcasted ephemeris : NAV_header,NAV_data generated from broadcasted ephemeris .n or .p
%% 		- precise orbits : replace NAV_header,NAV_data by sp3_header,sp3_data, and set options.nav to 'sp3'
%% - X0_base : [base receiver position] or [X;Y;Z;cdtr;cGGTO;cGPGL]
%%       If spp processing, set X0_base to NaN
%%	- epoch : computation epoch
%%	- options : structure containing processing options ( optional )
%% {
%%   X0 : approximated coordinates (column vector of 6 elements [X;Y;Z;cdtr;cGGTO;cGPGL]) if X0 is not defined, cut_off set to 0
%%        default : X0 = [0 ;0; 0; 0; 0; 0];
%%   const : constellations used ('G' = GPS, 'R' = Glonass, 'E' = Galileo, for multi-constellation concatenate chars)
%%        default : 'G'
%%   freq : type of used data [F1,F2,iono_free] 
%%        default : iono_free
%%         - 'F1' : use F1 frequency obs
%%         - 'F2' : use F2 frequency obs (or F5 for Galileo)
%%         - 'iono_free' : ionosphere-free combination
%%		   default : 'iono_free'
%%   iono : type of correction :
%%		   - 'klobuchar' : klobuchar modelization -> if  nav = 'brdc'
%%		   - 'none' : no correction (if freq = 'iono_free', iono set to 'none')
%%         default : 'none' 
%%   nav : type of orbits
%%		   - 'brdc' : broadcasted ephemeris
%%         - 'sp3' : precise orbits
%%         default : 'brdc'    
%%   cut_off : elevation cut off in rad
%%         default : 3*pi/180 rad
%%   DGNSS : DGNSS computation
%%         - 1 : yes
%%         - other : no DGNSS computation
%%         default : 0
%% }
%%
%%
%% Output
%% - G : cell tab (G1 if spp or G1 + G2 if DGNSS or phase (G1 = base, G2 = rover))
%%   nb_constellation * 32 lines, 8 columns
%%   column 1 : constellation
%%   column 2 : PRN
%%   column 3 : mjd (set to 0 if no satellite)
%%   column 4 : Code 1
%%   column 5 : Code 2
%%   column 6 : Iono-free combination with code observation
%%   column 7 : Phase 1
%%   column 8 : Phase 2
%%   column 9 : Iono-free combination with phase observation
%%   column 10 : Ephemeris of this satellite at mjd (if sp3, Eph.mjd set to 1)
%%   column 11 : structure containing intermediary results
%%      {
%%        PR =  21340831.1682333               : Code corrected by dte, TGD and dtrelat (m)
%%        PRC1 =  21001056.4770000             : Initial C1 code
%%        PRC2 =  21001061.9650000             : Initial C2 code
%%        PRC3 =  21001048.0735000             : Initial iono-free combination on code
%%        PhaseL1 =  21001066.9963885          : Initial L1 phase
%%        PhaseL2 =  21001079.5890304          : Initial L2 phase
%%        PhaseL3 =  21001047.7139056          : Initial iono-free combination on phase
%%        PR0 =  21331256.8538257              : Initial code (m)
%%        tr =  56442.9996527778               : reception time (mjd)
%%        travel_time =  0.0711534139188574    : travel time (s)
%%        te =  56442.9996519543               : emission time (mjd)
%%        dtrelat = -3.10256072047897e-09       : relativist correction (s)
%%        TGD = 0                              : TGD (s)
%%        dte =  3.19395777860622e-05          : satellite clock error (s)
%%        te_corr =  56442.9996519539          : corrected emission time (mjd) 
%%        Xs =  13482112.5310270               |
%%        Ys = -12886016.8650971               | : satellite position at te (m)
%%        Zs =  18861799.1754938               |
%%        Az =  4.91239303897012               : satellite azimuth (rad)
%%        Ele =  0.870337539193024             : satellite elevation (rad)
%%        h =  160.461182295345                : satellite elevation (m)
%%        dtropo =  3.04510732151988           : tropospheric correction (m)
%%        diono = 0                            : ionospheric correction (m)
%%        PR_corr =  21340828.1231260          : corrected pseudo-range (m)
%%        Xsr =  13482045.8232355              |
%%        Ysr = -12886086.6583699              | : satellite position at te, with rotation correction (m)
%%        Zsr =  18861799.1754938              |
%%      }
%% - result : structure containing final results
%%    {
%%      t =                                    : gpstime time structure
%%      {
%%        s1970 =  1369872000
%%        mjd =  56442
%%        jd =  2456442.50000000
%%        jd50 =  23160
%%        wk =  1742
%%        wsec =  345600
%%        yyyy =  2013
%%        yy =  13
%%        mon =  5
%%        dd =  30
%%        hh = 0
%%        min = 0
%%        sec = 0
%%        doy =  150
%%        wd =  4
%%        dsec = 0
%%        dy =  2013.40793976728
%%        GMST =  16.5112396504555
%%        EQEQ =  1.93995502140448e-04
%%        GAST =  16.5114336459577
%%      }
%%    
%%
%%
%%      PosSat =                               : satellite positions [Xs Ys Zs]
%%    
%%         11444831.05221041   19872840.57973917   13728071.85325433
%%         21785914.68839299  -14994024.60020332   -2170721.88255118
%%          3212056.70511915   16183670.97379925   20991792.25531062
%%         24942174.32843519   -6388801.36023886    5402212.53356499
%%    
%%      Dobs =                                 : corrected pseudo-ranges
%%    
%%         22802832.7032533
%%         24249925.7279791
%%         22795477.1178442
%%         21767854.7426687
%%    
%%      DobsC1 =                                 : C1
%%    
%%         22802832.7032533
%%         24249925.7279791
%%         22795477.1178442
%%         21767854.7426687
%%
%%      DobsC2 =                                 : C2
%%    
%%         22802832.7032533
%%         24249925.7279791
%%         22795477.1178442
%%         21767854.7426687
%%
%%      DobsC3 =                                 : C3
%%    
%%         22802832.7032533
%%         24249925.7279791
%%         22795477.1178442
%%         21767854.7426687
%%
%%      PobsL1 =                               : L1 obs
%%     
%%         21379867.5080223 
%%         22483230.1784891 
%%         23821262.8009240 
%%         20239064.0619077 
%%
%%      PobsL2 =                               : L2 obs
%%     
%%         21379867.5080223
%%         22483230.1784891
%%         23821262.8009240
%%         20239064.0619077
%%
%%      PobsL3 =                               : L3 obs
%%     
%%         21379867.5080223
%%         22483230.1784891
%%         23821262.8009240
%%         20239064.0619077
%%    
%%      sat_index =                            : satellite id ['constPRN',index in G] constPRN : Format A1I2
%%      
%%        {
%%          [1,1] = G14
%%          [2,1] = G20
%%          [3,1] = G22
%%          [4,1] = G32
%%          [1,2] =  14
%%          [2,2] =  20
%%          [3,2] =  22
%%          [4,2] =  32
%%        }
%%      
%%    
%%      ElevSat =                              : satellite elevation (rad)
%%    
%%         0.558558791057047
%%         0.244832779752556
%%         0.555984102770661
%%         0.693238540918065
%%
%%      AzSat =                                : satellite azimuth (rad)
%%    
%%         0.558558791057047
%%         0.244832779752556
%%         0.555984102770661
%%         0.693238540918065
%%
%%      Dtropo =                               : tropospheric delay correction (m)
%%          
%%         3.08312899894975 
%%         4.19965813067309  
%%         6.73638267466045   
%%         2.47054186459352  
%%    
%%      nb_GPS =  5                            : number of present GPS satellites at mjd
%%      nb_GLO = 0                             : number of present Glonass satellites at mjd
%%      nb_GAL = 0                             : number of present Galileo satellites at mjd
%%      nb_GPS_calc =  4                       : number of present GPS satellites at mjd used in computation (after cut_off)
%%      nb_GLO_calc = 0                        : number of present Glonass satellites at mjd used in computation (after cut_off)
%%      nb_GAL_calc = 0                        : number of present Galileo satellites at mjd used in computation (after cut_off)
%%      nb_sat = 4                             : number of total present satellites at mjd used in computation (after cut_off)
%%    }
%%
%% - result_LS_code_base / result_LS_code rover : LS code estimation result structure. See calc_LS_code() for details
%%     in case of spp : both are set to 'struct'
%%     in case of DGNSS : result_LS_code rover is set to 'struct'
%%
%% Information : use calc_preprocessing_spp/DGNSS/phase instead of this function
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%% Default parameters, output initialization and constants

%%% Default parameters
X0 = [0; 0; 0; 0; 0; 0]; %[X, Y, Z, cdtr, cGGTO, cGPGL]
iono = 'none'; % iono_free -> no more correction
nav = 'brdc';
const = 'G';
sat_num = 32;
cut_off = 3*pi/180;
freq = 'iono_free';
DGNSS = 0;

%%% output
G1=cell(0);
G2=cell(0);
result=cell(0);
result.PosSat = [];
result.Dobs = [];
result.DobsC1 = [];
result.DobsC2 = [];
result.DobsC3 = [];
result.PobsL1 = [];
result.PobsL2 = [];
result.PobsL3 = [];
result.sat_index = [];
result.AzSat = [];
result.ElevSat = [];
result.Dtropo = [];
result.nb_GPS = 0;
result.nb_GLO = 0;
result.nb_GAL = 0;
result.nb_GPS_calc = 0;
result.nb_GLO_calc = 0;
result.nb_GAL_calc = 0;
result.nb_sat = 0;
result.t.mjd = 0;

result_LS_code_base = struct;
result_LS_code_rover = struct;

% constants
c = 299792458.0;


%%% Options

% Get options if they are specified

if nargin==9

   [options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS] = set_options(options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS);
    
end
%%% End Options

% test of epoch
Nepoch = size(RNX_data_rover.G,3);
if(epoch>Nepoch || epoch<1)
	tool_print_info('Epoch > epoch_max or epoch < 1',3);
	return 
end

%%%%%% Type of computation : 
% 1 : spp
% 2 : DGNSS
% 3 : phase

if (~isfield(RNX_header_base,'VERSION') && DGNSS~=1) 
	tc = 1;
elseif (isfield(RNX_header_base,'VERSION') && DGNSS==1)
	tc = 2;
elseif (isfield(RNX_header_base,'VERSION') && DGNSS~=1)
	tc = 3;

else 
	tool_print_info('Calc_preprocessing : unknown computation type',3);
	return;
end

%%%%%% Observations loading

% nb_sat init
nb_sat = 0; % number of valid sat

nb_GPS = 0; % number of valid GPS sat
nb_GLO = 0; % number of valid GLONASS sat
nb_GAL = 0; % number of valid Galileo sat

nb_GPS_calc = 0; % number of valid GPS sat used in position computation
nb_GLO_calc = 0; % number of valid GLONASS sat used in position computation
nb_GAL_calc = 0; % number of valid Galileo sat used in position computation


print_epoch=1;
str_print_epoch = '';

%%%%%%  Cell tab filling

if (tc == 2 || tc == 3) % RNX header base is defined - 1 = rover, 2 = base
	[G2, index_sat, mjd_epoch, G1] = observation_loading(RNX_header_rover, RNX_data_rover, NAV_header, NAV_data, const, epoch, nav, RNX_header_base, RNX_data_base);
else % spp code peprocessing
	[G2, index_sat, mjd_epoch] = observation_loading(RNX_header_rover, RNX_data_rover, NAV_header, NAV_data, const, epoch, nav);
end	
  
% mjd in output : 0 if no obs
t = mjd_t(mjd_epoch);

if sum(index_sat>1)
	result.t=mjd_t(G2{index_sat(1),3}); % gpstime time structure
else
	result.t.mjd=0;
end  
  
%%%%%% Observation corrections and orbit computation
TO1 = NaN;
TO2 = NaN;
PRC = struct;

%%%%%% spp
if (tc==1)
	[G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = obs_correction_sat_pos(G2, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0, TO2);
end

%%%%%% DGNSS
if (tc==2)

	% base
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = obs_correction_sat_pos(G1, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0_base, TO1);
	
	% cut off
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = data_cut_off(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, cut_off);

	% LS estimation : PRC computation
	t1 = t.mjd*ones(size(PosSat1,1),1);
	options_LS.constraint=0.02; % 2 cm
	[result_LS1] = calc_LS_code(t1,PosSat1,Dobs1,ElevSat1,sat_index1,X0_base,options_LS);

	% output
	result_LS_code_base = result_LS1;
	
	% PRC
	if(result_LS_code_base.t~=0)
		PRC.corr = -result_LS1.V(:);
		PRC.index = result_LS1.sat_index(:,1);
	else
		PRC=struct;
	end
	
	% PRC correction on rover
	[G2, PosSat2, Dobs2, Dobs02, Pobs2 ElevSat2, AzSat2, Dtropo2, sat_index2] = obs_correction_sat_pos(G2, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0, TO2);
	
end

%%%%%% phase

if (tc == 3) 

	% phase preprocessing without dtr and time offsets
	% rover
	[G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = obs_correction_sat_pos(G2, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0, TO2);
	
	% base
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = obs_correction_sat_pos(G1, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0_base, TO1);

	% time offset computation and satellite position correction

	% rover
	% cut off
	[G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = data_cut_off(G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2, cut_off);

	t2 = t.mjd*ones(size(PosSat2,1),1);
	[result_LS2] = calc_LS_code(t2,PosSat2,Dobs2,ElevSat2,sat_index2,X0);
	result_LS_code_rover = result_LS2; % output 
	TO2 = [ result_LS2.cdtr/c ; result_LS2.cGGTO/c ; result_LS2.cGPGL/c ];
	[G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = obs_correction_sat_pos(G2, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0, TO2);

	% base
	% cut off
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = data_cut_off(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, cut_off);

	t1 = t.mjd*ones(size(PosSat1,1),1);
	options_LS_base.constraint=0.02; % 2 cm
	[result_LS1] = calc_LS_code(t1,PosSat1,Dobs1,ElevSat1,sat_index1,X0_base,options_LS_base);
	result_LS_code_base = result_LS1; % output 
	TO1 = [ result_LS1.cdtr/c ; result_LS1.cGGTO/c ; result_LS1.cGPGL/c ];
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = obs_correction_sat_pos(G1, index_sat, PRC, NAV_header, NAV_data, freq, nav,iono, X0_base, TO1);
end

% Cut off
[G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2 sat_index2] = data_cut_off(G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2, cut_off);

if (tc == 2 || tc == 3)
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1] = data_cut_off(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, cut_off);
end

% Check data (same satellites in base and rover stations)
if (tc == 2 || tc == 3)
	[G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = check_data(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2);
end

% sat num
for i = 1:size(sat_index2,1)
	cts = sat_index2{i};
	cts = cts(1);
	if strcmp(cts,'G')
		nb_GPS = nb_GPS + 1;
	elseif strcmp(cts,'R')
		nb_GLO = nb_GLO + 1;
	elseif strcmp(cts,'E')
		nb_GAL = nb_GAL + 1;
	end 
end
nb_sat = nb_GPS + nb_GLO + nb_GAL;

% print info
str_print_epoch = sprintf('PROCESSING EPOCH %04d-%02d-%02d %02d:%02d:%02d',t.yyyy,t.mon,t.dd,t.hh,t.min,t.sec);
tool_print_info(sprintf('%s   SAT : %2d GPS, %2d GLO, %2d GAL',str_print_epoch,nb_GPS,nb_GLO,nb_GAL),4);

% Output

% Sat number
for i = 1:size(sat_index2,1)
	if strcmp(G2{sat_index2{i,2},1},'E')
		nb_GAL_calc = nb_GAL_calc+1; % number of satellites used in final computation
	elseif strcmp(G2{sat_index2{i,2},1},'G')
		nb_GPS_calc = nb_GPS_calc+1; % number of satellites used in final computation
	elseif strcmp(G2{sat_index2{i,2},1},'R')
		nb_GLO_calc = nb_GLO_calc+1; % number of satellites used in final computation
	end
end

% rover if spp
if (tc == 1)
	result.PosSat = PosSat2;
	result.Dobs = Dobs2;
	result.DobsC1 = Dobs02(:,1);
	result.DobsC2 = Dobs02(:,2);
	result.DobsC3 = Dobs02(:,3);
	result.PobsL1 = Pobs2(:,1);
	result.PobsL2 = Pobs2(:,2);
	result.PobsL3 = Pobs2(:,3);
	result.AzSat = AzSat2;
	result.ElevSat = ElevSat2;
	result.Dtropo = Dtropo2;
	
% [base,rover] if DGNSS
elseif (tc == 2)
	%[base rover]
	result.PosSat = [PosSat1, PosSat2];
	result.Dobs = [Dobs1,Dobs2];
	result.DobsC1 = [Dobs01(:,1),Dobs02(:,1)];
	result.DobsC2 = [Dobs01(:,2),Dobs02(:,2)];
	result.DobsC3 = [Dobs01(:,3),Dobs02(:,3)];
	result.PobsL1 = [Pobs1(:,1),Pobs2(:,1)];
	result.PobsL2 = [Pobs1(:,2),Pobs2(:,2)];
	result.PobsL3 = [Pobs1(:,3),Pobs2(:,3)];
	result.AzSat = [AzSat1,AzSat2];
	result.ElevSat = [ElevSat1,ElevSat2];
	result.Dtropo = [Dtropo1,Dtropo2];
	

		
% [base,rover] if phase
elseif (tc == 3)
	%[base rover]
	result.PosSat = [PosSat1, PosSat2];
	result.Dobs = [Dobs1,Dobs2];
	result.DobsC1 = [Dobs01(:,1),Dobs02(:,1)];
	result.DobsC2 = [Dobs01(:,2),Dobs02(:,2)];
	result.DobsC3 = [Dobs01(:,3),Dobs02(:,3)];
	result.PobsL1 = [Pobs1(:,1),Pobs2(:,1)];
	result.PobsL2 = [Pobs1(:,2),Pobs2(:,2)];
	result.PobsL3 = [Pobs1(:,3),Pobs2(:,3)];
	result.AzSat = [AzSat1,AzSat2];
	result.ElevSat = [ElevSat1,ElevSat2];
	result.Dtropo = [Dtropo1,Dtropo2];
	
end

result.sat_index = sat_index2;

% sat number
% Before cut_off
result.nb_GPS = nb_GPS;
result.nb_GLO = nb_GLO;
result.nb_GAL = nb_GAL;

% After cut_off
result.nb_GPS_calc = nb_GPS_calc;
result.nb_GLO_calc = nb_GLO_calc;
result.nb_GAL_calc = nb_GAL_calc;

result.nb_sat = size(PosSat2,1);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







function [G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = check_data(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2)
%% function [G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2] = check_data(G1, PosSat1, Dobs1, Dobs01, Pobs1, ElevSat1, AzSat1, Dtropo1, sat_index1, G2, PosSat2, Dobs2, Dobs02, Pobs2, ElevSat2, AzSat2, Dtropo2, sat_index2)
%% Check than same satellites are visible from base and rover stations
%% Suppress obs of satellites visible from only one receiver
%%
%% Clement Fontaine 2013-12-18
%%
%% Input : 
%% - G1 :  cell tab containing satellite and PR informations (see calc_preprocessing for details)
%% - PosSat1 : matrix containing satellite positions (one line = one satellite)
%% - Dobs1 : vector containing corrected PR obs
%% - Dobs01 : vector containing initial PR obs [C1,C2,C3]
%% - Pobs1 : vector containing initial Phase obs [L1,L2,L3]
%% - ElevSat1 : vector containing satellite elevation
%% - AzSat1 : vector containing satellite azimuth
%% - Dtropo1 : vector containing tropospheric delay correction
%% - sat_index : cell containing index of satellites corresponding to observation in Pobs1 {constPRN, index}
%% 
%% idem for station 2
%%
%% Output : 
%% - idem than input, with size(PosSat1) == size(PosSat2)
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ind1 = cell2mat(sat_index1(:,2));
ind2 = cell2mat(sat_index2(:,2));

ind3 = sort([ind1;ind2]);
ind3 = ind3(sort(find(diff(ind3)==0)));

% Suppress in base
i = 1;
while i <=length(ind1)

	pos = find(ind3 == ind1(i));
	
	if length(pos) == 0
	
		PosSat1(i,:) = [];
		Dobs1(i,:) = [];
		Dobs01(i,:) = [];
		Pobs1(i,:) = [];
		ElevSat1(i,:) = [];
		AzSat1(i,:) = [];
		Dtropo1(i,:) = [];
		sat_index1(i,:) = [];
		ind1(i) = [];
	else
		i = i + 1;
	end
end


% Suppress in rover
i = 1;
while i <= length(ind2)

	pos = find(ind3 == ind2(i));

	if length(pos) == 0
		PosSat2(i,:) = [];
		Dobs2(i,:) = [];
		Dobs02(i,:) = [];
		Pobs2(i,:) = [];
		ElevSat2(i,:) = [];
		AzSat2(i,:) = [];
		Dtropo2(i,:) = [];
		sat_index2(i,:) = [];
		ind2(i) = [];

	else
		i = i + 1;
	end
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function [G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo sat_index] = data_cut_off(G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo, sat_index, cut_off)
%% function [G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo sat_index] = data_cut_off(G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo, sat_index, cut_off)
%% Suppress obs if satellite elevation < cut_off
%%
%% Clement Fontaine 2013-12-18
%%
%% Input : 
%% - G :  cell tab containing satellite and PR informations (see calc_preprocessing for details)
%% - PosSat : matrix containing satellite positions (one line = one satellite)
%% - Dobs : vector containing corrected PR obs
%% - Dobs0 : vector containing initial PR obs [C1,C2,C3]
%% - Pobs : vector containing initial Phase obs [L1,L2,L3]
%% - ElevSat : vector containing satellite elevation
%% - AzSat : vector containing satellite azimuth
%% - Dtropo : vector containing tropospheric delay correction
%% - sat_index : cell containing index of satellites corresponding to observation in Pobs1 {constPRN, index}
%% - cut_off : cut off value in rad
%%
%% Output : 
%% - idem than input, with only obs from satellites of which elevation < cut_off
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=1;
while i<=size(ElevSat,1)

	if(ElevSat(i)<cut_off)
	
		PosSat(i,:) = [];
		Dobs(i,:) = [];
		Dobs0(i,:) = [];
		Pobs(i,:) = [];
		ElevSat(i,:) = [];
		AzSat(i,:) = [];
		Dtropo(i,:) = [];
		sat_index(i,:) = [];
		
	else
		i = i + 1;	
	end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function [G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo, index_sat] = obs_correction_sat_pos(G, index_sat, PRC, NAV_header, NAV_data, freq, nav, iono, X0, TO)
%% function [G, PosSat, Dobs, Dobs0, Pobs, ElevSat, AzSat, Dtropo, index_sat] = obs_correction_sat_pos(G, index_sat, PRC, NAV_header, NAV_data, freq, nav, iono, X0, TO)
%% Corrected obs and satellite position computation
%%
%% Clement Fontaine 2013-12-18
%%
%% Input : 
%% - G : cell matrix containing loaded obs
%% - index_sat : index of loaded obs in G : G{index_sat(i),:} = obs of one satellite
%% - PRC : %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%TODO
%%  - NAV_header,NAV_data : nav data (rinex_n or sp3) generated by load_rinex_n (for broadcasted ephemeris) or load_sp3 (for precise orbits)
%% 		- broadcasted ephemeris : NAV_header,NAV_data generated from broadcasted ephemeris .n or .p
%% 		- precise orbits : replace NAV_header,NAV_data by sp3_header,sp3_data, and set options.nav to 'sp3'
%% - freq : used frequency ('F1', 'F2' or 'iono_free')
%% - nav : ephemeris type ('nav' for broadcasted eph, 'sp3' for precise eph)
%% - iono : ionospheric correction ('none' or 'klobuchar')
%% - X0 : approximated position of station
%% - TO : Time offset [dtr,GPGL,GGTO], in order to compute orbit in constellation time. Set to NaN if no correction of TO.
%%
%% Output : 
%% - G :  cell tab containing satellite and PR informations (see calc_preprocessing for details)
%% - PosSat : matrix containing satellite positions (one line = one satellite)
%% - Dobs : vector containing corrected PR obs
%% - Dobs0 : vector containing initial PR obs [C1,C2,C3]
%% - Pobs : vector containing initial Phase obs [L1,L2,L3]
%% - ElevSat : vector containing satellite elevation
%% - AzSat : vector containing satellite azimuth
%% - Dtropo : vector containing tropospheric delay correction
%% - sat_index : cell containing index of satellites corresponding to observation in Pobs1 {constPRN, index}. Format A1I2 : ex {'G02','R23'}
%%
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% constants
c = 299792458.0;
Rmin=6356000; % rayon min pour savoir si un point est sur terre
Rmax=6380000; % rayon max pour savoir si un point est sur terre

% limit values for PR
PRmin = 15000000;
PRmax = 35000000;

% meteo for tropospheric correction
T = 285; % Kelvin
P = 1013.15; % Hpa
H = 50; % percent

% frequencies
f1 = 1575.42;
f2 = 1227.60;
f5 = 1176.45;


nb_sat = size(index_sat,1);

% matrix initialization   
PosSat = zeros(nb_sat,3);
AzSat = zeros(nb_sat,1);
ElevSat = zeros(nb_sat,1);
Dobs = zeros(nb_sat,1);
Dobs0 = zeros(nb_sat,3);
Pobs = zeros(nb_sat,3);
Dtropo = zeros(nb_sat,1);

% index
index_sat_temp = zeros(nb_sat,1);


index_mat = 0;

% compute corrected PR and satellite position
for i=1:nb_sat
	
	numsat = index_sat(i);
		
	% frequency selection 
	if strcmp(freq,'iono_free')
		G{numsat,11}.PR = G{numsat,6};
	elseif strcmp(freq,'F1')
		G{numsat,11}.PR = G{numsat,4};
	elseif strcmp(freq,'F2')
		G{numsat,11}.PR = G{numsat,5};
	end
	
	% test iono_free + Glo + sp3
	if G{numsat,11}.PR==0 
		continue;
	end
	
	% raw observations
	G{numsat,11}.PRC1 = G{numsat,4};
	G{numsat,11}.PRC2 = G{numsat,5};
	G{numsat,11}.PRC3 = G{numsat,6};	
	G{numsat,11}.PhaseL1 = G{numsat,7};
	G{numsat,11}.PhaseL2 = G{numsat,8};
	G{numsat,11}.PhaseL3 = G{numsat,9};


		
	G{numsat,11}.PR0 = G{numsat,11}.PR;
	
	% Correction of PR in case of DGNSS
	if isfield(PRC,'corr')
	
		% if PR can be corrected
		[PR_correction] = corr_DGNSS(PRC,G{numsat,1},G{numsat,2});
		if PR_correction ~= 0
			G{numsat,11}.PR = G{numsat,11}.PR + PR_correction;
		else
			% no correction found
			continue;
		end
	
	end
	
	
	% TO already computed
    if(~isnan(TO)) % TO is defined
    
		% dtr
		G{numsat,11}.PR = G{numsat,11}.PR + c*TO(1);
	
		% GGTO
		if(strcmp(G{numsat,1},'E'))
			G{numsat,11}.PR = G{numsat,11}.PR + c*TO(2);
		end
	
		% GPGL
		if(strcmp(G{numsat,1},'R'))
			G{numsat,11}.PR = G{numsat,11}.PR + c*TO(3);
		end
		
		
    end
	
    % Reception date (modified julian date)
    G{numsat,11}.tr =  G{numsat,3};	
       

    % Emission date (modified julian date) 
    G{numsat,11}.travel_time = G{numsat,11}.PR / c; % travel_time
    G{numsat,11}.te = G{numsat,11}.tr - G{numsat,11}.travel_time / 86400 ; % te : emission time (mjd)

   
    % for Glonass : LEAP_SECONDS
    if(strcmp(G{numsat,1},'R') && strcmp(nav,'brdc'))
		G{numsat,11}.te = G{numsat,11}.te - NAV_header.LEAP_SECONDS/86400; % only if brdc
    end

    % relativist correction computation
    if strcmp(nav,'brdc')
		G{numsat,11}.dtrelat = corr_dtrelat_nav(G{numsat,10},G{numsat,11}.te);		
	else 
		G{numsat,11}.dtrelat = corr_dtrelat_sp3(NAV_data,G{numsat,1},G{numsat,2},G{numsat,11}.te,1e-3,9);
	end

    % TGD (Or BGD for Galileo)
    
    % On P1 : TGD
    % On P2 : (f1/f2)^2 * TGD
    % On P3 : 0
    
    if strcmp(nav,'brdc');
		% GPS
		if strcmp(G{numsat,1},'G')
		
			if strcmp(freq,'iono_free')
				G{numsat,11}.TGD = 0;
			elseif strcmp(freq,'F1')
				G{numsat,11}.TGD = -G{numsat,10}.TGD;
			elseif strcmp(freq,'F2')
				G{numsat,11}.TGD = -(f1/f2)^2*G{numsat,10}.TGD;
			end
	
		% Galileo % Combination between E1 and E5a
		elseif strcmp(G{numsat,1},'E')
			if strcmp(freq,'iono_free')
				G{numsat,11}.TGD = 0;
			elseif strcmp(freq,'F1')
				G{numsat,11}.TGD = -G{numsat,10}.BGDE5a;
			elseif strcmp(freq,'F2')
				G{numsat,11}.TGD = -(f1/f5)^2*G{numsat,10}.BGDE5a;
			end
					
		else
		
			G{numsat,11}.TGD = 0.0;
		end
		
	else
	
		G{numsat,11}.TGD = 0.0;
		
	end

    % Satellite clock drift
    if strcmp(nav,'brdc')
        G{numsat,11}.dte = corr_dte_nav(G{numsat,10},G{numsat,11}.te);
    else
		G{numsat,11}.dte = corr_dte_sp3(NAV_data,G{numsat,1},G{numsat,2},G{numsat,11}.te,9);
    end
    
    % Corrected emission time and corrected PR
    G{numsat,11}.te_corr = G{numsat,11}.te - G{numsat,11}.dte / 86400 - G{numsat,11}.TGD / 86400 - G{numsat,11}.dtrelat / 86400;    
    G{numsat,11}.PR = G{numsat,11}.PR +  c * G{numsat,11}.dte + c * G{numsat,11}.TGD + c * G{numsat,11}.dtrelat;
  
    % Satellite position computation
    if strcmp(nav,'brdc')
		[G{numsat,11}.Xs,G{numsat,11}.Ys,G{numsat,11}.Zs,dte] = orb_sat(G{numsat,10},G{numsat,1},G{numsat,2},G{numsat,11}.te_corr);
    else
		[G{numsat,11}.Xs,G{numsat,11}.Ys,G{numsat,11}.Zs,dte] = orb_sat(NAV_data,G{numsat,1},G{numsat,2},G{numsat,11}.te_corr,9);
    end
    
    
    
    
    
    % Az and Ele computation : for cut_off and tropospheric delay computation
    % Computation if initial coordinates are on the earth
    
    R = sqrt(sum(X0.^2));
    if (R>Rmin && R<Rmax) % on the earth ?
	    [Az,ele,h] = tool_az_ele_h(X0(1),X0(2),X0(3),G{numsat,11}.Xs,G{numsat,11}.Ys,G{numsat,11}.Zs);
	       
	    G{numsat,11}.Az = Az;
	    G{numsat,11}.Ele = ele;
	    G{numsat,11}.h = h;
	    
	    %~ % cut off
	    %~ if(ele < cut_off)
			 %~ continue
		%~ end
		
		% tropospheric delay - Saastamoinen model
		dtropo = corr_dtropo_saast(P,T,H,h,pi/2 - ele);
		
	else
		G{numsat,11}.Az = 0;
	    G{numsat,11}.Ele = 0;
	    G{numsat,11}.h = 0;
		dtropo = 11;	
	end
	
	G{numsat,11}.dtropo = dtropo;
	    
    % Ionospheric correction computation 

    if strcmp(freq,'iono_free') 
    
		G{numsat,11}.diono = 0; % no additional correction
	
	elseif (strcmp(freq,'F1') || strcmp(freq,'F2'))
	
		if (strcmp(iono,'klobuchar') && strcmp(nav,'brdc'))
			
			alpha_ion = NAV_header.GPSA;
			beta_ion = NAV_header.GPSB;
			G{numsat,11}.diono = corr_iono_klobuchar(X0(1),X0(2),X0(3),G{numsat,11}.Xs,G{numsat,11}.Ys,G{numsat,11}.Zs,alpha_ion,beta_ion,G{numsat,11}.te_corr,G{numsat,10},freq);
	    else
	    
			G{numsat,11}.diono = 0;
			
		end

	else
	
		G{numsat,11}.diono = 0;
		
	end
		
	% Ionospheric and tropospheric corrections 
    G{numsat,11}.PR_corr = G{numsat,11}.PR - G{numsat,11}.diono - G{numsat,11}.dtropo; 

    % time rotation during wave propagation
    alpha = -2 * pi * G{numsat,11}.PR_corr / 86400 / c;   
    [G{numsat,11}.Xsr,G{numsat,11}.Ysr,G{numsat,11}.Zsr] = tool_rotZ(G{numsat,11}.Xs,G{numsat,11}.Ys,G{numsat,11}.Zs,alpha);
	
	% matrix filling if non-0 position of satellite

	if (sqrt(G{numsat,11}.Xs^2+G{numsat,11}.Ys^2+G{numsat,11}.Zs^2)>PRmin && sqrt(G{numsat,11}.Xs^2+G{numsat,11}.Ys^2+G{numsat,11}.Zs^2)<PRmax)
		index_mat = index_mat+1;
		
		PosSat(index_mat,:) = [G{numsat,11}.Xsr,G{numsat,11}.Ysr,G{numsat,11}.Zsr];
		ElevSat(index_mat) = G{numsat,11}.Ele;
		AzSat(index_mat) = G{numsat,11}.Az;
		Dtropo(index_mat) = G{numsat,11}.dtropo;
		Dobs(index_mat) = G{numsat,11}.PR_corr;
		Dobs0(index_mat,:) = [G{numsat,11}.PRC1,G{numsat,11}.PRC2,G{numsat,11}.PRC3];
		Pobs(index_mat,:) = [G{numsat,11}.PhaseL1,G{numsat,11}.PhaseL2,G{numsat,11}.PhaseL3];
		index_sat_temp(index_mat) = index_sat(i);	
		
		
	end
		
end
% redim matrix
PosSat = PosSat(1:index_mat,:);
Dobs = Dobs(1:index_mat);
Dobs0 = Dobs0(1:index_mat,:);
Pobs = Pobs(1:index_mat,:);
Dtropo = Dtropo(1:index_mat,:);

AzSat = AzSat(1:index_mat); 
ElevSat = ElevSat(1:index_mat); 

index_sat = index_sat_temp(1:index_mat);

cell_index_sat = cell(length(index_sat),2);
for i = 1:length(index_sat)
	cell_index_sat{i,1} = sprintf('%s%02d',G{index_sat(i),1},G{index_sat(i),2});
	cell_index_sat{i,2} = index_sat(i);
end

index_sat = cell_index_sat;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






function [G1,index_sat, mjd, G2] = observation_loading(RNX_header1, RNX_data1, NAV_header, NAV_data, const, epoch, nav, RNX_header2, RNX_data2)
%% function [G1,index_sat, mjd, G2] = observation_loading(RNX_header1, RNX_data1, NAV_header, NAV_data, const, epoch, nav, RNX_header2, RNX_data2)
%% Load obs in a cell tab (one by obs)
%% If RNX 2 is defined, research sat in common between RNX1 and RNX2
%%
%% Clement Fotaine 2013-12-18
%%
%% Input : 
%% - RNX_header1, RNX_data1 : obs rinex data generated by function load_rinex_o
%%        Corresponds to rover station if RNX is defined
%% - NAV_header,NAV_data : nav data (rinex_n or sp3) generated by load_rinex_n (for broadcasted ephemeris) or load_sp3 (for precise orbits)
%% 		- broadcasted ephemeris : NAV_header,NAV_data generated from broadcasted ephemeris .n or .p
%% 		- precise orbits : replace NAV_header,NAV_data by sp3_header,sp3_data, and set options.nav to 'sp3'
%% - const : const to load : 'G' for GPS, 'R' for Glonass and 'E' for Galileo. Several constellations : 'GRE' or 'GE' or 'GR'
%% - epoch : epoch to which obs are going to be loaded
%% - nav : ephemeris type ('nav' for broadcasted eph, 'sp3' for precise eph)
%% - RNX_header2, RNX_data2 : obs rinex data generated by function load_rinex_o (optional)
%%        Corresponds to base station if defined

%% Output : 
%% - G1 : cell matrix containing loaded obs of RNX 1
%% - index_sat : index of loaded obs in G : G{index_sat(i),:} = obs of one satellite
%% - mjd : date of obs in Modified Julian Day
%% - G2 : cell matrix containing loaded obs of RNX 2 (optional, if RNX is defined)
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mjd = 0;

% limit values for PR
PRmin = 15000000;
PRmax = 35000000;

% if differential : RNX1 = rover and RNX2 = base

% cell tab init
sat_num = 32*length(const);
G1 = cell(sat_num,11);  % data
index_sat = zeros(sat_num,1);  % index of valid sat
nb_sat = 0; % number of valid sat

if nargin == 9 % differencial positionning -> two observation structures

	G2 = cell(sat_num,11);

end

for i=1:length(const) % n constellations

	current_const = const(i);
	
	for j = 1:32 % 32 satellites by constellation
			
		current_num = 32*(i-1)+j;
		
		G1{current_num,1} = current_const; %    constellation
		G1{current_num,2} = j; %                PRN
		 
		G1{current_num,3} = 0; %                mjd (set to 0 if no obs)         
		G1{current_num,4} = 0; %                C1
		G1{current_num,5} = 0; %                P2
		G1{current_num,6} = 0; %                P3
		G1{current_num,7} = 0; %                L1
		G1{current_num,8} = 0; %                L2
		G1{current_num,9} = 0; %                L3
		G1{current_num,10} = 0; %               Eph
		G1{current_num,11} = struct; %          Structure containing corrected PR and satellite positions
		
		
		% Get observations in rnx1		
		Obs1 = get_obs(RNX_header1,RNX_data1,current_const,j,epoch);
		
		if isfield(Obs1,'mjd') % there is this satellite at this epoch
		
			mjd = Obs1.mjd;
			
			% Get ephemeris of this satellite
			if(strcmp(nav,'brdc'))

				[Eph]=get_ephemeris(NAV_header,NAV_data,current_const,j,mjd);

			elseif(strcmp(nav,'sp3')) % NAV_data = sp3_data

				% test : approximated position in order to check if sp3 contain the satellite
				Eph = [];
				[X_test,Y_test,Z_test,dte_test] = orb_sat(NAV_data,current_const,j,mjd,1);              
				if(~sum(X_test)==0)
					Eph.mjd=1;
				end
								
			end
			
			if isfield(Eph,'mjd')
								
				% Obs : depends on the constellation
				
				[PR1,PR2,Phase1,Phase2] = get_select_obs(Obs1,Eph);

				
				% missing information
				if(PR1==0 || PR2==0 || Phase1==0 || Phase2==0) 
					continue
				end
				
				% test PR
				if(PR1>PRmax || PR1<PRmin || PR2>PRmax || PR2<PRmin) 
					continue
				end
				
				% differential GNSS -> seek corresponding epoch in base rinex
				if nargin == 9
				
					G2{current_num,1} = current_const; %    constellation
					G2{current_num,2} = j; %                PRN
		 
					G2{current_num,3} = 0; %                mjd (set to 0 if no obs)         
					G2{current_num,4} = 0; %                C1
					G2{current_num,5} = 0; %                P2
					G2{current_num,6} = 0; %                P3
					G2{current_num,7} = 0; %                L1
					G2{current_num,8} = 0; %                L2
					G2{current_num,9} = 0; %                L3
					G2{current_num,10} = 0; %               Eph
					G2{current_num,11} = struct; %          Structure containing corrected PR and satellite position
				
					epoch_base = get_epoch_from_mjd(RNX_header2,mjd);
					
					Obs2 = get_obs(RNX_header2,RNX_data2,current_const,j,epoch_base);
					
					[PR1_2,PR2_2,Phase1_2,Phase2_2] = get_select_obs(Obs2,Eph);
					
					% missing information
					if(PR1_2==0 || PR2_2==0 || Phase1_2==0 || Phase2_2==0) 
						continue
					end
				
					% test PR
					if(PR1_2>PRmax || PR1_2<PRmin || PR2_2>PRmax || PR2_2<PRmin) 
						continue
					end
				
				end		
				
				% fill G
				
				% Rover
				G1{current_num,3} = mjd; % mjd of observation
				
				G1{current_num,4} = PR1; % C1 (for GPS)
				G1{current_num,5} = PR2; % P2
				Eph.const = G1{current_num,1};
				Eph.PRN = G1{current_num,2};
				G1{current_num,6} = corr_iono_free(PR1,PR2,Eph); %%%%%%%%%%%%%%%%%%%%%% Glonass -> no iono free if sp3...
				G1{current_num,7} = Phase1; % C1 (for GPS)
				G1{current_num,8} = Phase2; % P2	
				G1{current_num,9} = corr_iono_free(Phase1,Phase2,Eph); %%%%%%%%%%%%%%%%%%%%%% Glonass -> no iono free if sp3...

				G1{current_num,10} = Eph; % Ephemeris					
				% Base
				if nargin == 9
					G2{current_num,3} = mjd; % mjd of observation
					
					G2{current_num,4} = PR1_2; % C1 (for GPS)
					G2{current_num,5} = PR2_2; % P2
					G2{current_num,6} = corr_iono_free(PR1_2,PR2_2,Eph); %%%%%%%%%%%%%%%%%%%%%% Glonass -> no iono free if sp3...
					G2{current_num,7} = Phase1_2; % C1 (for GPS)
					G2{current_num,8} = Phase2_2; % P2	
					G2{current_num,9} = corr_iono_free(Phase1_2,Phase2_2,Eph); %%%%%%%%%%%%%%%%%%%%%% Glonass -> no iono free if sp3...
	
					G2{current_num,10} = Eph; % Ephemeris	
				end
					
				nb_sat = nb_sat + 1;
				index_sat(nb_sat) = current_num;
			
			else 		
				% can't find Ephemeris of this satellite			
			end
					
		else
			% No satellite at this epoch
		end	
		 
	end
	
end

index_sat = index_sat(1:nb_sat);


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function [PR1,PR2,Phase1,Phase2] = get_select_obs(Obs,Eph)
%% function [PR1,PR2,Phase1,Phase2] = select_obs(Obs,Eph)
%% Load obs from obs structure for a given satellite (in Eph)
%%
%% Clement Fontaine 2013-12-18
%%
%% Input : 
%% - Obs : obs structure (from get_obs())
%% - Eph : eph structure (from get_ephemeris(), or just with fields 'const' ans 'PRN' defined)
%%
%% Output : 
%% - PR1 : first pseudo-range obs
%% - PR2 : second pseudo-range obs
%% - Phase1 : first phase obs
%% - Phase2 : second phase obs
%%
%% See code for details about observation type 
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = 299792458.0;

% output
PR1 = 0;
PR2 = 0;
Phase1 = 0;
Phase2 = 0;

% Obs : depends on the constellation
if isfield(Obs,'constellation')
	const = Obs.constellation;
else
	return;
end

if strcmp(const,'G') % GPS
	% RNX v2.11 -> C1, P2, L1, L2
	% RNX v3.02 -> C1C, C2W, L1C, L2W
	f1 = 1575.42;
	f2 = 1227.60;
	
	PR1 = Obs.C1;
	PR2 = Obs.C2;
	
	Phase1 = Obs.L1*c/(f1*1e6);
	Phase2 = Obs.L2*c/(f2*1e6);
	

elseif strcmp(const,'R') % GLONASS
	% RNX v2.11 -> C1, P2, L1, L2
	% RNX v3.02 -> C1C, C2C, L1C, L2C
	
	PR1 = Obs.C1;
	PR2 = Obs.C2;
	
	if(nargin == 2)
		if isfield(Eph,'freq_num')
			f1 = 1602.00 + Eph.freq_num*0.5625;
			f2 = 1246.00 + Eph.freq_num*0.4375;
			Phase1 = Obs.L1*c/(f1*1e6);
			Phase2 = Obs.L2*c/(f2*1e6);
		else
			tool_print_info('Glonass constellation with sp3 file : freq_num is missing, Phase1 = Phase2 = 1',2);
			f1 = 1;
			f2 = 1;
			Phase1 = 1;
			Phase2 = 1;
		end
	else
		tool_print_info('Glonass constellation with sp3 file : freq_num is missing, Phase1 = Phase2 = 1',2);
		f1 = 1;
		f2 = 1;
		Phase1 = 1;
		Phase2 = 1;
	end

elseif(strcmp(const,'E')) % Galileo
	% RNX v2.11 -> C1, C5, L1, L5
	% RNX v3.02 -> C1C or C1X ,C5Q or C5X, L1C or L1X, L5Q or L5X 
	f1 = 1575.42;
	f2 = 1176.45;
		
	PR1 = Obs.C1;
	PR2 = Obs.C2;	
	
	Phase1 = Obs.L1*c/(f1*1e6);
	Phase2 = Obs.L2*c/(f2*1e6);	

end


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




function [options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS] = set_options(options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS)
%% function [options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS] = set_options(options,X0,freq,cut_off,iono,nav,const,sat_num,DGNSS)
%% Set options
%% 
%% Clement Fontaine 2013-12-18
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	 
if isfield(options,'freq')
	if (strcmp(options.freq,'F1') || strcmp(options.freq,'F2') || strcmp(options.freq,'iono_free')) 
		freq = options.freq;  
	else    
		options.freq = 'iono_free';
	end 
else    
	options.freq = 'iono_free';
end 

if isfield(options,'cut_off')
	cut_off = options.cut_off; % input in rad
else    
	options.cut_off = cut_off;
end

% keep following code after cut off
if isfield(options,'X0')

	if options.X0(1) == 0
	
		tool_print_info('No approximates values, cut_off set to 0',2);
		options.X0 = X0;
		cut_off = 0;
		
	else
	
		X0 = options.X0(:); % in column	

	end	

else
	tool_print_info('No approximates values, cut_off set to 0',2);
	options.X0 = X0;
	cut_off = 0;
end  


if isfield(options,'iono') 
	if (strcmp(options.iono,'klobuchar') || strcmp(options.iono,'none')) 
		iono = options.iono; 
	else 
		options.iono = iono;
	end
else 
	options.iono = iono;
end    

if isfield(options,'nav') 
	if (strcmp(options.nav,'sp3'))
		nav = 'sp3'; 
	else 
		options.nav = nav;
	end  
else 
	options.nav = nav;
end  

if isfield(options,'DGNSS') 
	if (options.DGNSS==1)
		DGNSS = 1; 
	else 
		options.DGNSS = DGNSS;
	end  
else 
	options.DGNSS = DGNSS;
end 

if isfield(options,'const')
	sat_num = 0;
	const = '';
	if regexp(options.const,'G')>0
		const = strcat(const,'G');
		sat_num = sat_num + 32;
	end
	if regexp(options.const,'R')>0
		const = strcat(const,'R');
		sat_num = sat_num + 32;
	end	
	if regexp(options.const,'E')>0
		const = strcat(const,'E');
		sat_num = sat_num + 32;
	end      
else 
	options.const = const;
end  

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
